## Replication Code
## Accountability for the Local Economy at All Levels of Government in United States Elections
## Justin de Benedictis-Kessner Christopher Warshaw
## This code file creates all maps

library(foreign)
library(arm)
library(lfe)
library(readr)
library(stargazer)
library(tidyverse)

lag.new <- function(x, n = 1L, along_with){
	index <- match(along_with - n, along_with, incomparable = NA)
	out <- x[index]
	attributes(out) <- attributes(x)
	out
}

## load data
load(file="econ_counties_cities_analysis.Rdata")

data_analysis$pres_dem_control2 <- NA
data_analysis$pres_dem_control2[data_analysis$pres_dem_control==0] <- -1
data_analysis$pres_dem_control2[data_analysis$pres_dem_control==1] <- 1

## merge NA counties from VA with cities according to BEA merge patterns from here https://apps.bea.gov/regional/docs/statelist.cfm
## Have to bring in VA FIPS codes for joined city/county areas:
va_counties <- read_csv("VA_cities_counties_link2.csv")
va_counties <- va_counties %>%
	mutate(fips = as.character(fips)) %>%
	spread(
		key = "county",
		value="fips") %>%
	rename(city_fips = `0`,
				 county_fips = `1`)

va_cities <- filter(va_counties,!is.na(city_fips))
va_counties <- filter(va_counties,!is.na(county_fips))
va_counties$city_fips <- va_cities$city_fips[match(va_counties$bea_fips,va_cities$bea_fips)]


# View(filter(data_analysis,state_abb=="VA"))
# View(filter(data_analysis, year==2004 & is.na(wages_perworker_real_delta)))
# missing: some in VA; Shannon County, SD; Maui, HI

# setting groups BEA county codes back to their real FIPS codes:
data_analysis$fips[which(data_analysis$fips %in% va_counties$bea_fips)] <- va_counties$county_fips[match(data_analysis$fips[which(data_analysis$fips %in% va_counties$bea_fips)],va_counties$bea_fips)]
# View(filter(data_analysis,state_abb=="VA" & year==2008))


## two-way FE +LDV
width=.25
counties_pop2010 <- subset(data_analysis,year==2010) %>%
	as_tibble() %>%
	dplyr::select(fips,population) %>%
	rename(population_2010 = population) %>%
	filter(!is.na(population_2010))
data_analysis <- left_join(data_analysis,counties_pop2010,by="fips")

data_analysis$random <- rnorm(nrow(data_analysis))
data_analysis$wages_perworker_real_delta_account <- data_analysis$wages_perworker_real_delta*data_analysis$pres_dem_control2

## Democratic presidents
#reg2.graph.outcome2<-(lm(pres_demshare~factor(state_year), data=data_analysis[!is.na(data_analysis$pres_demshare)& !is.na(data_analysis$wages_perworker_real_delta),]))
#reg2.graph.treatment2<-(lm(wages_perworker_real_delta_account~factor(state_year), weights= data_analysis $pres_total_votes_scaled[!is.na(data_analysis$pres_total_votes_scaled) & !is.na(data_analysis$pres_demshare)], data=data_analysis[!is.na(data_analysis$pres_total_votes_scaled)& !is.na(data_analysis$pres_demshare),]))

reg2.graph.outcome2 <- (felm(pres_demshare_delta ~ random | fips + state_year, 
														 data=subset(data_analysis,!is.na(pres_demshare) & 
														 							# population_2010 > 20000 &
														 							!is.na(wages_perworker_real_delta))
														 ))
# reg2.graph.treatment2 <- (felm(wages_perworker_real_delta_account ~ random | fips + state_year,
# 															 data=subset(data_analysis,!is.na(pres_demshare) & 
# 															 							# population_2010 > 20000 & 
# 															 							!is.na(pres_demshare_delta))
# 															 ))

reg2.graph.treatment2 <- (felm(wages_perworker_real_delta ~ random | fips + state_year,
															 data=subset(data_analysis,!is.na(pres_demshare) & 
															 							# population_2010 > 20000 & 
															 							!is.na(pres_demshare_delta))
))

# extract residualized DV and IV:
predictions.outcome2 <- resid(reg2.graph.outcome2)
predictions.treatment2 <- resid(reg2.graph.treatment2)

data.predictions2 <- subset(data_analysis,!is.na(wages_perworker_real_delta) &
															!is.na(pres_demshare_delta) &!is.na(pres_demshare) &
														!is.na(fips) 	&!is.na(state_year)&!is.na(random) #&
															# !is.na(population_2010) &
															# population_2010 > 20000
														)
colnames(predictions.outcome2) <- "predictions.outcome"
colnames(predictions.treatment2) <- "predictions.treatment"
data.predictions2 <- cbind(data.predictions2, predictions.outcome2, predictions.treatment2)

reg21 <- felm(pres_demshare ~ wages_perworker_real_delta_account | fips + state_year,
							data=subset(data_analysis,
													# population_2010 > 20000 &
													!is.na(pres_demshare) 
														)
)
reg22 <- (lm(predictions.outcome ~ predictions.treatment, data=data.predictions2))
summary(reg21)
summary(reg22)

StateCodes <- read_csv("StateCodes.csv")
StateCodes <- dplyr::select(StateCodes, Name,POAbrv)
# data.predictions2 <- dplyr::select(data.predictions2, county_name.x, state_abb,predictions.treatment, year,predictions.outcome,county_newspaper) 
data.predictions2 <- dplyr::select(data.predictions2, fips, county_name.x, state_abb,year,wages_perworker_real_delta,predictions.treatment, predictions.outcome) 
data.predictions2 <- merge(data.predictions2, StateCodes, by.x="state_abb", by.y="POAbrv")

# data.predictions2$predictions.treatment2 <- NA
# data.predictions2$predictions.treatment2[data.predictions2$predictions.treatment< -.33] <- 0
# data.predictions2$predictions.treatment2[data.predictions2$predictions.treatment> -.33 & data.predictions2$predictions.treatment<.33] <- 1
# data.predictions2$predictions.treatment2[data.predictions2$predictions.treatment>.33] <- 2

data.predictions2$predictions.treatment2 <- NA
data.predictions2$predictions.treatment2[data.predictions2$predictions.treatment<0] <- 0
data.predictions2$predictions.treatment2[data.predictions2$predictions.treatment>0] <- 1

sum(!is.na(data.predictions2$predictions.treatment2[which(data.predictions2$year==1992)]))
sum(!is.na(data.predictions2$predictions.treatment2[data.predictions2$year==2004]))
sum(!is.na(data.predictions2$predictions.treatment2[data.predictions2$year==2008]))
sum(!is.na(data.predictions2$predictions.treatment2[data.predictions2$year==2012]))
# data.predictions2 %>% filter(year==2012 & state_abb=="MA") %>% View()
# data_analysis %>% filter(year==2012 & state_abb=="MA") %>% View()

summary(reg22)
reg22 <- (lm(predictions.outcome~factor(predictions.treatment2), data=data.predictions2))




#data.predictions2$fips[data.predictions2$fips=="12086"]<-  "12025"


### Maps
library(ggplot2)
library(ggmap)
library(sf)
library(tigris)
library(tidycensus)
options(tigris_class = "sf",tigris_use_cache = T)

# states <- map_data("state")
# counties <- map_data("county")
# counties <- tigris::counties(cb = T,resolution = "20m",shift_geo=T)
counties <- get_estimates(geography = "county",
													variables = "BIRTHS",
															 geometry = TRUE,
															 shift_geo = TRUE) # this shifts Hawaii and Alaska over

counties$merge_fips <- counties$GEOID
counties <- mutate(counties,
									 state_fips = str_sub(GEOID,1,2))

# counties <- map_data("city")
# LA <- map_data("state", region="Maryland")
# View(filter(data.predictions2, year==2008 & state_abb=="VA"))
cities <- read_csv(file="big_cities.csv")

# counties$subregion <- gsub(" ","", counties$subregion)
# data.predictions2$county_name.x <- gsub(" ","", data.predictions2$county_name.x)

# counties2 <- geo_join(spatial_data = counties, data_frame =  data.predictions2, by_sp="GEOID", by_df="fips", how="left")


counties2 <- left_join(counties,data.predictions2,by=c("GEOID" = "fips"))
# sanity checks:
# View(filter(counties2, state_abb=="VA"))
# View(filter(counties, !(GEOID %in% counties2$GEOID[which(counties2$year==2004)])))


# counties2 <- arrange(counties2, order)
# counties2 <- subset(counties2, year==2008 )
# (counties_base <- ggplot(counties) + 
# 	geom_polygon(data=counties,col="black",fill="white",aes(x=long,y=lat))
# 	)

ggplot(counties) + 
	geom_sf(data=counties,col="grey30",lwd=0.5,inherit.aes=F)


counties_continuous <- ggplot(counties) + 
	geom_sf(data=counties,lwd=0.15,col="grey30",fill="white") + 
	geom_sf(data=subset(counties2,year==1992),aes(fill = cut_interval(predictions.treatment,6)),col="grey30",lwd=0,inherit.aes=F,show.legend = T) +
						scale_fill_brewer(palette="RdYlGn",name="Change in wages\nrelative to state-year avg.") + 
						# theme(legend.key.width = unit(0.75, "cm"),legend.key.height = unit(2, "cm"),
						# 			axis.line = NULL,panel.border = NULL,panel.grid = NULL) +
						# theme_bw()
						theme(axis.line=element_blank(),axis.text.x=element_blank(),
									axis.text.y=element_blank(),axis.ticks=element_blank(),
									axis.title.x=element_blank(),
									axis.title.y=element_blank(),
									panel.background=element_blank(),panel.border=element_blank(),panel.grid.major = element_line(colour = 'transparent'),
									panel.grid.minor=element_blank(),plot.background=element_blank())

#pdf("Figures/counties_continuous_wages_1992.pdf", height=3.75, width=7.5)
#counties_continuous
#dev.off()

counties_continuous <- ggplot(counties) + 
	geom_sf(data=counties,lwd=0.15,col="grey30",fill="white") + 
	geom_sf(data=subset(counties2,year==2004),aes(fill = cut_interval(predictions.treatment,6)),col="grey30",lwd=0,inherit.aes=F,show.legend = T) +
	scale_fill_brewer(palette="RdYlGn",name="Change in wages\nrelative to state-year avg.") + 
	# theme(legend.key.width = unit(0.75, "cm"),legend.key.height = unit(2, "cm"),
	# 			axis.line = NULL,panel.border = NULL,panel.grid = NULL) +
	# theme_bw()
	theme(axis.line=element_blank(),axis.text.x=element_blank(),
				axis.text.y=element_blank(),axis.ticks=element_blank(),
				axis.title.x=element_blank(),
				axis.title.y=element_blank(),
				panel.background=element_blank(),panel.border=element_blank(),panel.grid.major = element_line(colour = 'transparent'),
				panel.grid.minor=element_blank(),plot.background=element_blank())

#pdf("Figures/counties_continuous_wages_2004.pdf", height=3.75, width=7.5)
#counties_continuous
#dev.off()

counties_continuous <- ggplot(counties) + 
	geom_sf(data=counties,lwd=0.15,col="grey30",fill="white") + 
	geom_sf(data=subset(counties2,year==2008),aes(fill = cut_interval(predictions.treatment,6)),col="grey30",lwd=0,inherit.aes=F,show.legend = T) +
	scale_fill_brewer(palette="RdYlGn",name="Change in wages\nrelative to state-year avg.") + 
	# theme(legend.key.width = unit(0.75, "cm"),legend.key.height = unit(2, "cm"),
	# 			axis.line = NULL,panel.border = NULL,panel.grid = NULL) +
	# theme_bw()
	theme(axis.line=element_blank(),axis.text.x=element_blank(),
				axis.text.y=element_blank(),axis.ticks=element_blank(),
				axis.title.x=element_blank(),
				axis.title.y=element_blank(),
				panel.background=element_blank(),panel.border=element_blank(),panel.grid.major = element_line(colour = 'transparent'),
				panel.grid.minor=element_blank(),plot.background=element_blank())

#pdf("Figures/counties_continuous_wages_2008.pdf", height=3.75, width=7.5)
#counties_continuous
#dev.off()


counties_continuous <- ggplot(counties) + 
	geom_sf(data=counties,lwd=0.15,col="grey30",fill="white") + 
	geom_sf(data=subset(counties2,year==2012),aes(fill = cut_interval(predictions.treatment,6)),col="grey30",lwd=0,inherit.aes=F,show.legend = T) +
	scale_fill_brewer(palette="RdYlGn",name="Change in wages\nrelative to state-year avg.") + 
	# theme(legend.key.width = unit(0.75, "cm"),legend.key.height = unit(2, "cm"),
	# 			axis.line = NULL,panel.border = NULL,panel.grid = NULL) +
	# theme_bw()
	theme(axis.line=element_blank(),axis.text.x=element_blank(),
				axis.text.y=element_blank(),axis.ticks=element_blank(),
				axis.title.x=element_blank(),
				axis.title.y=element_blank(),
				panel.background=element_blank(),panel.border=element_blank(),panel.grid.major = element_line(colour = 'transparent'),
				panel.grid.minor=element_blank(),plot.background=element_blank())

#pdf("Figures/counties_continuous_wages_2012.pdf", height=3.75, width=7.5)
#counties_continuous
#dev.off()

counties_binary <- ggplot(counties) + 
	geom_sf(data=counties,col="grey30",fill="white",lwd=0.15) + 
	geom_sf(data=subset(counties2,year==1976),aes(fill = predictions.treatment2),col="grey30",lwd=0,inherit.aes=F,show.legend = T) +
	scale_fill_gradient(name="Change in wages\nrelative to state-year avg.",
											low="#d73027", high="#1a9850", guide='legend', breaks=c(0,1),
											labels=c("shrinking","growing")) +
	# theme(legend.key.width = unit(0.75, "cm"),legend.key.height = unit(2, "cm"),
	# 			axis.line = NULL,panel.border = NULL,panel.grid = NULL) +
	# theme_bw()
	theme(axis.line=element_blank(),axis.text.x=element_blank(),
				axis.text.y=element_blank(),axis.ticks=element_blank(),
				axis.title.x=element_blank(),
				axis.title.y=element_blank(),
				panel.background=element_blank(),panel.border=element_blank(),panel.grid.major = element_line(colour = 'transparent'),
				panel.grid.minor=element_blank(),plot.background=element_blank())


pdf("Figures/counties_binary_wages_1976.pdf", height=3.75, width=7.5)
counties_binary
dev.off()

counties_binary <- ggplot(counties) + 
	geom_sf(data=counties,col="grey30",fill="white",lwd=0.15) + 
	geom_sf(data=subset(counties2,year==1992),aes(fill = predictions.treatment2),col="grey30",lwd=0,inherit.aes=F,show.legend = T) +
	scale_fill_gradient(name="Change in wages\nrelative to state-year avg.",
											low="#d73027", high="#1a9850", guide='legend', breaks=c(0,1),
											labels=c("shrinking","growing")) +
	# theme(legend.key.width = unit(0.75, "cm"),legend.key.height = unit(2, "cm"),
	# 			axis.line = NULL,panel.border = NULL,panel.grid = NULL) +
	# theme_bw()
	theme(axis.line=element_blank(),axis.text.x=element_blank(),
				axis.text.y=element_blank(),axis.ticks=element_blank(),
				axis.title.x=element_blank(),
				axis.title.y=element_blank(),
				panel.background=element_blank(),panel.border=element_blank(),panel.grid.major = element_line(colour = 'transparent'),
				panel.grid.minor=element_blank(),plot.background=element_blank())


pdf("Figures/counties_binary_wages_1992.pdf", height=3.75, width=7.5)
counties_binary
dev.off()

counties_binary <- ggplot(counties) + 
	geom_sf(data=counties,col="grey30",fill="white",lwd=0.15) + 
	geom_sf(data=subset(counties2,year==2004),aes(fill = predictions.treatment2),col="grey30",lwd=0,inherit.aes=F,show.legend = T) +
	scale_fill_gradient(name="Change in wages\nrelative to state-year avg.",
											low="#d73027", high="#1a9850", guide='legend', breaks=c(0,1),
											labels=c("shrinking","growing")) +
	# theme(legend.key.width = unit(0.75, "cm"),legend.key.height = unit(2, "cm"),
	# 			axis.line = NULL,panel.border = NULL,panel.grid = NULL) +
	# theme_bw()
	theme(axis.line=element_blank(),axis.text.x=element_blank(),
				axis.text.y=element_blank(),axis.ticks=element_blank(),
				axis.title.x=element_blank(),
				axis.title.y=element_blank(),
				panel.background=element_blank(),panel.border=element_blank(),panel.grid.major = element_line(colour = 'transparent'),
				panel.grid.minor=element_blank(),plot.background=element_blank())


pdf("Figures/counties_binary_wages_2004.pdf", height=3.75, width=7.5)
counties_binary
dev.off()

counties_binary <- ggplot(counties) + 
	geom_sf(data=counties,col="grey30",fill="white",lwd=0.15) + 
	geom_sf(data=subset(counties2,year==2008),aes(fill = predictions.treatment2),col="grey30",lwd=0,inherit.aes=F,show.legend = T) +
	scale_fill_gradient(name="Change in wages\nrelative to state-year avg.",
											low="#d73027", high="#1a9850", guide='legend', breaks=c(0,1),
											labels=c("shrinking","growing")) +
	# theme(legend.key.width = unit(0.75, "cm"),legend.key.height = unit(2, "cm"),
	# 			axis.line = NULL,panel.border = NULL,panel.grid = NULL) +
	# theme_bw()
	theme(axis.line=element_blank(),axis.text.x=element_blank(),
				axis.text.y=element_blank(),axis.ticks=element_blank(),
				axis.title.x=element_blank(),
				axis.title.y=element_blank(),
				panel.background=element_blank(),panel.border=element_blank(),panel.grid.major = element_line(colour = 'transparent'),
				panel.grid.minor=element_blank(),plot.background=element_blank())


pdf("Figures/counties_binary_wages_2008.pdf", height=3.75, width=7.5)
counties_binary
dev.off()

counties_binary <- ggplot(counties) + 
	geom_sf(data=counties,col="grey30",fill="white",lwd=0.15) + 
	geom_sf(data=subset(counties2,year==2012),aes(fill = predictions.treatment2),col="grey30",lwd=0,inherit.aes=F,show.legend = T) +
	scale_fill_gradient(name="Change in wages\nrelative to state-year avg.",
											low="#d73027", high="#1a9850", guide='legend', breaks=c(0,1),
											labels=c("shrinking","growing")) +
	# theme(legend.key.width = unit(0.75, "cm"),legend.key.height = unit(2, "cm"),
	# 			axis.line = NULL,panel.border = NULL,panel.grid = NULL) +
	# theme_bw()
	theme(axis.line=element_blank(),axis.text.x=element_blank(),
				axis.text.y=element_blank(),axis.ticks=element_blank(),
				axis.title.x=element_blank(),
				axis.title.y=element_blank(),
				panel.background=element_blank(),panel.border=element_blank(),panel.grid.major = element_line(colour = 'transparent'),
				panel.grid.minor=element_blank(),plot.background=element_blank())


pdf("Figures/counties_binary_wages_2012.pdf", height=3.75, width=7.5)
counties_binary
dev.off()

### Greyscale versions:
counties_binary <- ggplot(counties) + 
	geom_sf(data=counties,col="grey30",fill="white",lwd=0.15) + 
	geom_sf(data=subset(counties2,year==1992),aes(fill = predictions.treatment2),col="grey30",lwd=0,inherit.aes=F,show.legend = T) +
	scale_fill_gradient(name="Change in wages\nrelative to state-year avg.",
											low="lightgrey", high="grey23", guide='legend', breaks=c(0,1),
											labels=c("shrinking","growing")) +
	# theme(legend.key.width = unit(0.75, "cm"),legend.key.height = unit(2, "cm"),
	# 			axis.line = NULL,panel.border = NULL,panel.grid = NULL) +
	# theme_bw()
	theme(axis.line=element_blank(),axis.text.x=element_blank(),
				axis.text.y=element_blank(),axis.ticks=element_blank(),
				axis.title.x=element_blank(),
				axis.title.y=element_blank(),
				panel.background=element_blank(),panel.border=element_blank(),panel.grid.major = element_line(colour = 'transparent'),
				panel.grid.minor=element_blank(),plot.background=element_blank())


pdf("Figures/counties_binary_wages_greyscale_1992.pdf", height=3.75, width=7.5)
counties_binary
dev.off()

counties_binary <- ggplot(counties) + 
	geom_sf(data=counties,col="grey30",fill="white",lwd=0.15) + 
	geom_sf(data=subset(counties2,year==2004),aes(fill = predictions.treatment2),col="grey30",lwd=0,inherit.aes=F,show.legend = T) +
	scale_fill_gradient(name="Change in wages\nrelative to state-year avg.",
											low="lightgrey", high="grey23", guide='legend', breaks=c(0,1),
											labels=c("shrinking","growing")) +
	# theme(legend.key.width = unit(0.75, "cm"),legend.key.height = unit(2, "cm"),
	# 			axis.line = NULL,panel.border = NULL,panel.grid = NULL) +
	# theme_bw()
	theme(axis.line=element_blank(),axis.text.x=element_blank(),
				axis.text.y=element_blank(),axis.ticks=element_blank(),
				axis.title.x=element_blank(),
				axis.title.y=element_blank(),
				panel.background=element_blank(),panel.border=element_blank(),panel.grid.major = element_line(colour = 'transparent'),
				panel.grid.minor=element_blank(),plot.background=element_blank())


pdf("Figures/counties_binary_wages_greyscale_2004.pdf", height=3.75, width=7.5)
counties_binary
dev.off()

counties_binary <- ggplot(counties) + 
	geom_sf(data=counties,col="grey30",fill="white",lwd=0.15) + 
	geom_sf(data=subset(counties2,year==2008),aes(fill = predictions.treatment2),col="grey30",lwd=0,inherit.aes=F,show.legend = T) +
	scale_fill_gradient(name="Change in wages\nrelative to state-year avg.",
											low="lightgrey", high="grey23", guide='legend', breaks=c(0,1),
											labels=c("shrinking","growing")) +
	# theme(legend.key.width = unit(0.75, "cm"),legend.key.height = unit(2, "cm"),
	# 			axis.line = NULL,panel.border = NULL,panel.grid = NULL) +
	# theme_bw()
	theme(axis.line=element_blank(),axis.text.x=element_blank(),
				axis.text.y=element_blank(),axis.ticks=element_blank(),
				axis.title.x=element_blank(),
				axis.title.y=element_blank(),
				panel.background=element_blank(),panel.border=element_blank(),panel.grid.major = element_line(colour = 'transparent'),
				panel.grid.minor=element_blank(),plot.background=element_blank())


pdf("Figures/counties_binary_wages_greyscale_2008.pdf", height=3.75, width=7.5)
counties_binary
dev.off()

counties_binary <- ggplot(counties) + 
	geom_sf(data=counties,col="grey30",fill="white",lwd=0.15) + 
	geom_sf(data=subset(counties2,year==2012),aes(fill = predictions.treatment2),col="grey30",lwd=0,inherit.aes=F,show.legend = T) +
	scale_fill_gradient(name="Change in wages\nrelative to state-year avg.",
											low="lightgrey", high="grey23", guide='legend', breaks=c(0,1),
											labels=c("shrinking","growing")) +
	# theme(legend.key.width = unit(0.75, "cm"),legend.key.height = unit(2, "cm"),
	# 			axis.line = NULL,panel.border = NULL,panel.grid = NULL) +
	# theme_bw()
	theme(axis.line=element_blank(),axis.text.x=element_blank(),
				axis.text.y=element_blank(),axis.ticks=element_blank(),
				axis.title.x=element_blank(),
				axis.title.y=element_blank(),
				panel.background=element_blank(),panel.border=element_blank(),panel.grid.major = element_line(colour = 'transparent'),
				panel.grid.minor=element_blank(),plot.background=element_blank())


pdf("Figures/counties_binary_wages_greyscale_2012.pdf", height=3.75, width=7.5)
counties_binary
dev.off()


